##############################################################################
# File-Name: 02_analysis_primary_outcomes.r
# author: Tiago Ventura
# Purpose: Analysis primary outcomes
# Machine: macOS Ventura 13.5.1
# R version 4.3.1 
##############################################################################

# Packages ----------------------------------------------------------------
pacman::p_load(tidyverse, scales,  knitr, srvyr, extrafont, 
               rebus, broom, tidyr, lubridate, here, ggdist, ggtext, ggalt, janitor, 
               broom, tidyr,  stargazer, janitor, estimatr, list, marginaleffects, xtable)


# Source Auxiliary Code ---------------------------------------------------
source(here("code","_utils.r"))


# Data -----
d <- readRDS(here("data","all_processed_data.rds"))


# Figure 3 -------------------------------------------------

## Dependent variable: sum of recall to false rumors
depvar = "false_news_exp"
## See models at _utils.r
false_false_e_itt <- func_ols_base(depvar, data=d)
false_false_e_ittc <- func_ols_covariates(depvar, data=d)
false_false_e_cace <- func_cace(depvar, data=d)



## Dependent variable: sum of exposure to true news
depvar = "true_news_e"

# models
true_true_e_itt <- func_ols_base(depvar, d)
true_true_e_ittc <- func_ols_covariates(depvar,d)
true_true_e_cace <- func_cace(depvar, data=d)

# plot all together
output= list(false_false_e_itt, false_false_e_ittc, false_false_e_cace,
             true_true_e_itt, true_true_e_ittc, true_true_e_cace)

models= list("ITT", "Cov-ITT", "CACE", "ITT", "Cov-ITT", "CACE")

dep_var_label= unlist(list(rep("False Rumors Exposure", 3), rep("True News Exposure", 3)))

# plot
plot_model(output, models, dep_var_label, 
           title="A",
           caption = "Standardized Treatment Effects \n on Exposure to False and True Headlines")

save_function("fig3.png")
save_function("fig3.tiff")


# Figure 4  ---------------------------------------------------------------

## Dependent variable: misinformation beleifs
depvar = "false_false_sum"

#models
false_false_itt <- func_ols_base(depvar, data=d)
false_false_ittc <- func_ols_covariates(depvar, data=d)
false_false_cace <- func_cace(depvar, data=d)

## Dependent variable: Mainstream News Knowledge: sum of say true news are true
depvar = "true_true_sum"

# models
true_true_itt <- func_ols_base(depvar, data=d)
true_true_ittc <- func_ols_covariates(depvar,d)
true_true_cace <- func_cace(depvar, data=d)

# plot all together
output= list(false_false_itt, false_false_ittc, false_false_cace,
             true_true_itt, true_true_ittc, true_true_cace)

models= list("ITT", "Cov-ITT", "CACE", "ITT", "Cov-ITT", "CACE")
dep_var_label= unlist(list(rep("False Rumors Accuracy", 3), rep("True News Accuracy", 3)))

# plot
plot_model(output, models, dep_var_label,
           title="B",
           caption = "Standardized Treatment Effects \n on Belief Accuracy on False and True Headlines")

save_function("fig4.png")
save_function("fig4.tiff")

# Fig 5 --------------------------------------------------------------------

## data wrangling for moderator

# convert whatsapp usage for politics to numneric
d$whats_pol_n <- as.numeric(fct_rev(d$w1_q_whatsapp_purposes_politics))

# consider don't know responses as NA
d$whats_pol_n <- ifelse(d$whats_pol_n==1, NA,d$whats_pol_n-2)

# separate by terciles
quantile(d$whats_pol_n, na.rm = TRUE)

d <- d %>% 
    mutate(whats_pol_heavy=ifelse(whats_pol_n>3, "First Tercile", 
                                  ifelse(whats_pol_n==0, "Third Tercile", NA)), 
           whats_pol_heavy=fct_rev(whats_pol_heavy))


# model for False news belief + covariates
model <- lm_robust(false_false_sum ~ exp*whats_pol_heavy + 
                       
                       #demographics
                       w1_q_gender +  
                       q_race +  
                       w1_q_education + 
                       w1_q_age+
                       income_num +

                       
                       # trust  
                       as.numeric(w1_q_trust_government) +
                       as.numeric(w1_q_trust_congress) +
                       as.numeric(w1_q_trust_supreme_court) +
                       as.numeric(w1_q_trust_electoral_authorities) +
                       as.numeric(w1_q_trust_globo) +
                       as.numeric(w1_q_trust_news_channels) +
                       w1_q_legitimacy + 
                       
                       # polarization    
                       w1_affective_pol_bolsonaro +
                       w1_affective_pol_lula +

                       # ideology  
                       w1_ideology_you +  
                       w1_q_politics_num +   
                       
                       # social media usage
                       w1_q_whatsapp_num +
                       as.numeric(w1_q_fake_news) +
                       as.numeric(w1_q_whatsapp_purposes_news) +
                       as.numeric(w1_q_whatsapp_purposes_pay_bills) +
                       as.numeric(w1_q_whatsapp_purposes_family) +
                       as.numeric(w1_q_whatsapp_purposes_groups), data=d)




# get p-value
stars <- car::linearHypothesis(model, c("exptreatment:whats_pol_heavyFirst Tercile +exptreatment =  exptreatment"))
pval = stars$`Pr(>Chisq)`[2]

# get marginal effects with binary moderator
false <- interaction_plot_binary(model, effect="exptreatment", 
                                 moderator="whats_pol_heavyFirst Tercile", 
                                 interaction="exptreatment:whats_pol_heavyFirst Tercile",
                                 conf=.95) 


# graph
false <- false %>%
    mutate(moderator_f=fct_rev(ifelse(moderator==0, "Low use of WhatsApp \n to share/receive political \n information", 
                                      "High use of WhatsApp \n to share/receive political \n information")), 
           lb95=lower_bound,
           up95=upper_bound, 
           up90=delta_1 + 1.64*std.error,
           lb90=delta_1 - 1.64*std.error)

# Plot Marginal Effects Interactive
ggplot(false, 
       aes(y=delta_1, 
           color=moderator_f, 
           x=moderator_f)) +
    geom_errorbar(aes(ymin=lb90, 
                      ymax=up90), 
                  size=2, width=0, alpha=.8, 
                  position=position_dodge(width = .8), 
                  color=colorwes) +
    geom_errorbar(aes(ymin=lb95, 
                      ymax=up95),  
                  size=1, width=.1, alpha=.5, 
                  position=position_dodge(width = .8), 
                  color=colorwes) +
    geom_point(size=14, shape=21, fill="white",
               position=position_dodge(width = .6), 
               stroke=2, 
               color=colorwes, alpha=1) +
    labs(title="", 
         y="Marginal Effect", 
         x="", 
         caption="") +
    geom_segment(aes(x=1.2, xend=1.8, y=0.30, yend=0.30), 
                 size=.8, 
                 color="black") +
    geom_segment(aes(x=1.2, xend=1.2, y=0.30, yend=0.28), 
                 size=.8, 
                 color="black") +
    geom_segment(aes(x=1.8, xend=1.8, y=0.30, yend=0.28), 
                 size=.8, 
                 color="black") +
    geom_text(aes(x=1.5, y=.35, 
                  label=paste0("p-value=", round(pval, digits=3)), family=""),
              size=4, family=my_font, color = "black")  +
    geom_text(aes(label=round(delta_1, digits=2)), 
              fontface = "italic", 
              fill="white", 
              color = "black",
              size=4, 
              position = position_dodge(.6))

save_function("fig5b.png")
save_function("fig5b.tiff")

# Estimate the model with continuous moderator 
model <- lm_robust(false_false_sum ~ exp*whats_pol_n + 
                       
                       #demographics
                       w1_q_gender +  
                       q_race +  
                       w1_q_education + 
                       w1_q_age+  
                       income_num +
                       
                       # trust  
                       as.numeric(w1_q_trust_government) +
                       as.numeric(w1_q_trust_congress) +
                       as.numeric(w1_q_trust_supreme_court) +
                       as.numeric(w1_q_trust_electoral_authorities) +
                       as.numeric(w1_q_trust_globo) +
                       as.numeric(w1_q_trust_news_channels) +
                       w1_q_legitimacy + 
                       
                       # polarization    
                       w1_affective_pol_bolsonaro +
                       w1_affective_pol_lula +
                       
                       # ideology  
                       w1_ideology_you +  
                       w1_q_politics_num +   
                       
                       # social media usage
                       w1_q_whatsapp_num +
                       as.numeric(w1_q_fake_news) +
                       as.numeric(w1_q_whatsapp_purposes_news) +
                       as.numeric(w1_q_whatsapp_purposes_pay_bills) +
                       as.numeric(w1_q_whatsapp_purposes_family) +
                       as.numeric(w1_q_whatsapp_purposes_groups), data=d)


# get marginal effects with continous
false <- interaction_plot_continuous_ggplot(model, effect="exptreatment", 
                                            moderator="whats_pol_n", 
                                            interaction="exptreatment:whats_pol_n",
                                            conf=.95, num_points = 10) 

ggplot(false, 
       aes(y=delta_1, ymax=upper_bound, ymin=lower_bound, x=moderator)) +
    geom_errorbar(size=1, 
                  width=.2, 
                  alpha=.6, 
                  color=colorwes) +
    geom_point(size=4, shape=21, fill="white",
               position=position_dodge(width = .6), 
               stroke=2, 
               color=colorwes, 
               alpha=1) +
    geom_ribbon(alpha=.1,
                color=NA,
                fill=colorwes,
                linetype="dashed")  +
    geom_hline(aes(yintercept=0), linetype="dashed", color=colorwes) +
    labs(title="", 
         y="Marginal Effect", 
         x=" \n Use of WhatsApp to share/receive \n political information") +
    scale_x_continuous(breaks = 0:5,
                       labels = c("Every few weeks", "", "", "", "", "At least 10 times \n a day")) +
    theme(axis.title.x = element_text(hjust=0.5))

save_function("fig5a.png")
save_function("fig5a.tiff")

# Fig 6 -------------------------------------------------------------------

## Dependent variable: Polarization ------------------------------------------------------------

# variable names
list_politems_items <- list("pol_all_zcore", 
                            "distance_candidates",
                            "affective_pol_diff_all", 
                            "social_pol_agg", 
                            "mean_abs_policy_polarization")
# Models
list_pol_itt <- map(list_politems_items, ~ func_ols_base(.x, data=d))
list_pol_ittc <- map(list_politems_items, ~ func_ols_covariates(.x, data=d))
list_pol_cace <- map(list_politems_items, ~ func_cace(.x, data=d))

# plot all together
output = flatten(list(list_pol_itt, list_pol_ittc, list_pol_cace))
models= c(rep("ITT", 5), rep("Cov-ITT", 5), rep("CACE", 5))
dep_var_label= rep(c("Polarization \n Index", 
                     "Ideological \n Polarization",
                     "Affective \n Polarization", 
                     "Social \n Polarization", 
                     "Issue \n Polarization"), 3)

# plots
s <- plot_model(output,
                models, 
                dep_var_label, 
                title="", 
                subtitle= "", 
                caption="Standardized Treatment Effects on Polarization") +
    theme(axis.text.x  = element_text(size=10), 
          strip.text = element_text(family = my_font, 
                                    color = "#22211d",
                                    size = 10, face="bold")
    )

# change the names of the outcomes
s$data <- s$data %>%
    mutate(dv=fct_relevel(dv, "Polarization \n Index", 
                          "Ideological \n Polarization",
                          "Affective \n Polarization", 
                          "Social \n Polarization", 
                          "Issue \n Polarization")) 


## DV: subjective well-being ---------------------------------------------------

# variable names
list_wb_items <- list("swb_zcore", # zscore index
                      "q_well_being_happy_numeric"  ,
                      "q_well_being_depressed_numeric",
                      "q_well_being_anxious_numeric",
                      "q_well_being_isolated_numeric", 
                      "q_well_being_satisfied_numeric" )

# run the models
list_swb_itt <- map(list_wb_items, ~ func_ols_base(.x, data=d))
list_swb_ittc <- map(list_wb_items, ~ func_ols_covariates(.x, data=d))
list_swb_cace <- map(list_wb_items, ~ func_cace(.x, data=d))

# plot all together
output = flatten(list(list_swb_itt, list_swb_ittc, list_swb_cace))
models= c(rep("ITT", 6), rep("Cov-ITT", 6), rep("CACE", 6))
dep_var_label= rep(c("Index", 
                     "Happy",
                     "Depressed", 
                     "Anxious", 
                     "Isolated", 
                     "Satisfied"), 3)


# get only the composite z-score
list_wb_items <- list("swb_zcore")
list_itt <- map(list_wb_items, ~ func_ols_base(.x, data=d))
list_ittc <- map(list_wb_items, ~ func_ols_covariates(.x, data=d))
list_cace <- map(list_wb_items, ~ func_cace(.x, data=d))

# plot all together
output = flatten(list(list_itt, list_ittc, list_cace))
models= c(rep("ITT", 1), rep("Cov-ITT", 1), rep("CACE", 1))
dep_var_label= rep(c(""), 3)

# Combine Polarization and SW ---------------------------------------------
list_bind_itt <- c(list_pol_itt, list(list_swb_itt[[1]]))
list_bind_ittc <- c(list_pol_ittc, list(list_swb_ittc[[1]]))
list_bind_cace <- c(list_pol_cace, list(list_swb_cace[[1]]))

# plot all together
output = flatten(list(list_bind_itt, list_bind_ittc,list_bind_cace))
models= c(rep("ITT", 6), rep("Cov-ITT", 6), rep("CACE", 6))
dep_var_label= rep(c("Polarization Index", 
                     "Ideological \n Polarization",
                     "Affective Polarization", 
                     "Social Polarization", 
                     "Issue Polarization", 
                     "Subjective Well-Being Index"), 3)

# only index
index <- plot_model(output,
                    models, 
                    dep_var_label, 
                    title="", 
                    caption ="") 

# all polarization
pol<- plot_model(output,
                 models, 
                 dep_var_label, 
                 title="", 
                 caption ="") 

# select only the two indexes
index$data <- index$data %>%
    filter(str_detect(dv, "Index"))

index +
    labs(caption="Standardized Treatment Effects on Polarization \n and Subjective Well-Being Index", 
         title="")


save_function("fig6.png")
save_function("fig6.tiff")




# Table 17 - Multiple hypothesis Adjustment ---------------------------------------------


# Selecting all p-values into a vector
P_Values_Hyp <- c(false_false_e_itt$p.value, 
                  true_true_itt$p.value,
                  true_true_itt$p.value, 
                  list_pol_itt[[1]]$p.value, 
                  list_swb_itt[[1]]$p.value)

# four hypothesis
P_Values_Hyp <- round(P_Values_Hyp,4)

#Using the Bwnjamini and Hochberg (1995) method report new p-values.
FDR_Adj_P <- p.adjust(P_Values_Hyp,method='fdr', n=length(P_Values_Hyp))
FDR_Adj_P <- as.numeric(FDR_Adj_P)


Hyp_list <- c('(H1)',
              '(H2a)',
              '(H2b)',
              '(H3)',
              '(H4)')

#Convert from numeric and character:
P_Values_Hyp <- as.character(P_Values_Hyp)
FDR_Adj_P <- as.character(FDR_Adj_P)

# build a matrix
mat_P <- matrix(c(P_Values_Hyp,
                  FDR_Adj_P),ncol=5,byrow=T)

# define names of the matrix
rownames(mat_P) <- c('Unadjusted P-Value',
                     'P-Value (FDR Adjusted)')
colnames(mat_P) <- Hyp_list

# save table
print(xtable(mat_P,
             caption='Unadjusted and Adjusted P-Values Testing Each Hypothesis'),
      caption.placement = 'top',
      file="output/sm_tab17.tex")


# Table 18 - MH:Polarization outcomes ------------------------------------------------------------

# get all p-values for polarization outcomes
P_Values_Hyp <- map(list_pol_itt, "p.value")


# adjustment method
FDR_Adj_P <- p.adjust(P_Values_Hyp,method='fdr', n=length(P_Values_Hyp))

# build matrix
mat_P <- matrix(c(P_Values_Hyp,
                  FDR_Adj_P),ncol=5,byrow=T)

rownames(mat_P) <- c('Unadjusted P-Value',
                     'P-Value (FDR Adjusted)')


colnames(mat_P) <- c("Polarization Index", "False Polarization",
                     "Affective Polarization", "Social Polarization", "Policy Polarization")


# save table
print(xtable(mat_P,
             caption='Unadjusted and Adjusted P-Values for Polarization Outcomes'),
      caption.placement = 'top',
      file="output/sm_tab18.tex")


#  Table 19 - MH:Subjective well-being ------------------------------------------------------------

# get all p-values for swb outcomes
P_Values_Hyp <- map(list_swb_itt, "p.value")

# adjustment
FDR_Adj_P <- p.adjust(P_Values_Hyp,method='fdr', n=length(P_Values_Hyp))


# build a matrix
mat_P <- matrix(c(P_Values_Hyp,
                  FDR_Adj_P),ncol=6,byrow=T)

rownames(mat_P) <- c('Unadjusted P-Value',
                     'P-Value (FDR Adjusted)')

colnames(mat_P) <- c("Index", 
                     "Happy"  ,
                     "Depressed",
                     "Anxious",
                     "Isolated", 
                     "Satisfied" )

# save table
print(xtable(mat_P,
             caption='Unadjusted and Adjusted P-Values for Subjective Wellbeing Outcomes'),
      caption.placement = 'top',
      file="output/sm_tab19.tex")



